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ABSTRACT 

The mass spectrum of stellar-mass black holes (BHs) is highly uncertain. Dynamical 
mass measurements are available only for few (~ 10) BHs in X-ray binaries, while the¬ 
oretical models strongly depend on the hydrodynamics of supernova (SN) explosions 
and on the evolution of massive stars. In this paper, we present and discuss the mass 
spectrum of compact remnants that we obtained with SEVN, a new public population- 
synthesis code, which couples the PARSEC stellar evolution tracks with up-to-date 
recipes for SN explosion (depending on the Carbon-Oxygen mass of the progenitor, 
on the compactness of the stellar core at pre-SN stage, and on a recent two-parameter 
criterion based on the dimensionless entropy per nucleon at pre-SN stage). SEVN can 
be used both as a stand-alone code and in combination with direct-summation N-body 
codes (Starlab, HiGPUs). The PARSEC stellar evolution tracks currently imple¬ 
mented in SEVN predict significantly larger values of the Carbon-Oxygen core mass 
with respect to previous models. For most of the SN recipes we adopt, this implies 
substantially larger BH masses at low metallicity 2x 10 3 ), than other population- 
synthesis codes. The maximum BH mass found with SEVN is ~ 25, 60 and 130 Mq 
at metallicity Z = 2 x 10 -2 , 2 x 10 -3 and 2 x 10 -4 , respectively. Mass loss by stellar 
winds plays a major role in determining the mass of BHs for very massive stars (^90 
Mq), while the remnant mass spectrum depends mostly on the adopted SN recipe for 
lower progenitor masses. We discuss the implications of our results for the transition 
between NS and BH mass, and for the expected number of massive BHs (with mass 
> 25 Mg) as a function of metallicity. 

Key words: black hole physics - methods: numerical - stars: black holes - stars: 
evolution - stars: mass-loss - stars: neutron 


1 INTRODUCTION 

Compact remnants are the final stage of the evolution of 
massive stars, and power a plethora of important astrophys- 
ic.al processes: they are the engine of the X-ray binaries we 
observe in the nearby Universe, and may be powerful sources 
of gravitational waves (e.g. Phinney 1991). Furthermore, the 
merger of two neutron stars (NSs) and/or that of a stellar 
black hole (BH) with a NS are expected to lead to one of 
the most energetic transient phenomena in the Universe: 
the short gamma-ray bursts (e.g. Paczynski 1991). Finally, 
X-ray binaries powered by BHs and/or NSs are the key to 
explain some of the most luminous point-like non-nuclear 
X-ray sources (the ultraluminous X-ray sources, e.g. Mapelli 
et al. 2010; Mapelli & Zampieri 2014 and references therein), 
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and are an important source of feedback, in both the nearby 
and the early Universe (e.g. Justham & Schawinski 2012, 
and references therein). 

Despite their importance for astrophysics, the details 
of the formation of BHs and NSs (and especially the link 
with their progenitor stars) are matter of debate. From the 
observational point of view, the confirmed BHs are only a 
few tens (see table 2 of Ozel et al. 2010, for one of the most 
updated compilations). These are located in X-ray binaries, 
mostly in the Milky Way (MW), and an accurate dynamical 
mass estimate has been derived only for a fraction of them 
(~ 10). Most of the derived BH masses are in the range 5 ^ 
tubh/ Mq ^ 10. In the MW, the most massive BHs in X-ray 
binaries do not significantly exceed m bh ~ 15 Mq, whereas 
a few BHs in nearby galaxies might have higher masses: M33 
X-7 (m B H = 15.65 ± 1.45 M 0 , Orosz et al. 2007), IC-10 X- 
1 (tobh ~ 23 — 34M 0 , Prestwich et al. 2007; Silverman & 
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Filippenko 2008), NGC 300 X-l (?tibh > 10 M@, Crowther 
et al. 2007, 2010). Interestingly, these three massive BHs 
are in regions with relatively low metallicity. A metallicity 
Z ~ 0.004 is estimated for the dwarf irregular galaxy IC-10 
(Garnett 1990). The metallicity of M33 in proximity of X-7 
is Z ~ 0.008, and that of NGC300 in proximity of X-l is 
Z ~ 0.006 (Pilyugin et al. 2004). 

The statistics is significantly larger for NSs: currently, 
there are dynamical mass measurements for 61 NSs (17, 
11, 30, and 3 of them are in X-ray binaries, NS-NS bina¬ 
ries, NS-white dwarf binaries and NS-main sequence bina¬ 
ries, respectively, http://stellarcollapse.org/nsmasses, 
Lattimer & Prakash 2005; Lattimer 2012). 

The link between the progenitor star and the compact 
remnant is still poorly constrained for both BHs and NSs: 
observations of core-collapse supernovae (SNe) indicate a 
deficit of massive (> 20 Mg) progenitor stars (Smartt 2009; 
Horiuchi et al. 2011; Jennings et al. 2012, 2014; Gerke et al. 

2014) , which possibly suggests that the most massive stars 
undergo no or faint SNe. 

From a theoretical perspective, the formation and the 
mass spectrum of BHs and NSs strongly depend on two fun¬ 
damental processes: (i) the hydrodynamics of SNe; (ii) mass 
loss by stellar winds in massive stars (during and especially 
after the main sequence, MS). 

(i) The physics of SN explosions is extremely complex, 
and the hydrodynamical codes that investigate the explosion 
mechanisms are computationally challenging (see e.g. Fryer 
1999; Fryer & Kalogera 2001; Heger & Woosley 2002; Heger 
et al. 2003; Fryer 2006; O’Connor & Ott 2011; Fryer et al. 
2012; Janka 2012; Ugliano et al. 2012; Burrows 2013; Pejcha 
& Prieto 2015; Ertl et al. 2015). In particular, the link be¬ 
tween the late evolutionary stages of a massive star and the 
SN products is still matter of debate. Several authors (e.g. 
Bethe 1990; Janka et al. 2007; Burrows 2013; Janka 2012) 
investigate for which structural properties of the progenitor 
star a SN can fail, leading to the direct collapse of the star 
to a BH. Even if the SN occurs, how much matter can fall 
back and be accreted onto the proto-compact remnant is 
very uncertain. 

(ii) For massive progenitors (zero-age MS mass 
Mzams ^ 30 Mq) the details of stellar evolution are very im¬ 
portant for the SN outcome and for the final remnant mass. 
In fact, the final mass Ma n of the progenitor star (i.e. the 
mass of a star immediately before the collapse) is governed 
by the amount of mass loss by stellar winds (e.g. Mapelli 
et al. 2009; Belczynski et al. 2010; Fryer et al. 2012; Mapelli 
et al. 2013). The rate of mass loss by stellar winds on the 
MS increases with the metallicity of the star as M oc Z a , 
where a ~ 0.5 — 0.9, depending on the model (e.g., Ku- 
dritzki et al. 1987; Leitherer et al. 1992; Kudritzki & Puls 
2000; Vink et al. 2001a; Kudritzki 2002). The behaviour of 
evolved massive stars, such as luminous blue variable stars 
(LBVs) and Wolf-Rayet stars (WRs), is also expected to 
depend on metallicity, but with larger uncertainties (e.g., 
Vink & de Koter 2005; Meynet & Maeder 2005; Bressan 
et al. 2012; Tang et al. 2014). 

Both the models of SN explosion (e.g. Fryer et al. 2012; 
Janka 2012; Burrows 2013; Pejcha & Prieto 2015; Ertl et al. 

2015) and the theory of massive star evolution (e.g. Bres¬ 
san et al. 2012; Tang et al. 2014) were deeply revised in the 
last few years. For these reasons, population synthesis codes 


that aim at studying the demographics of compact remnants 
must account for up-to-date models for both SN explosions 
and stellar evolution. Here we present SEVN (acronym for 
‘Stellar EVolution N-body’), a new population synthesis tool 
that couples PARSEC evolutionary tracks for stellar evolu¬ 
tion (Bressan et al. 2012; Chen et al. 2014; Tang et al. 2014) 
with up-to-date models for SN explosion (Fryer et al. 2012; 
Janka 2012; Ertl et al. 2015), and that can be easily merged 
with several N-body codes. The new PARSEC evolution¬ 
ary tracks consider the most recent updates for mass loss 
by stellar winds and other input physics. In this paper, we 
present and discuss the mass spectrum of BHs and NSs that 
we obtain from SEVN, with particular attention to the de¬ 
pendence of the remnant mass on metallicity. 

Furthermore, SEVN is extremely versatile, because it 
relies upon a set of tables extracted from stellar evolu¬ 
tion tracks: if we are interested in comparing different stel¬ 
lar evolution models, we can do it quickly and easily, by 
changing tables. The new tool is publicly available 1 . SEVN 
is specifically designed to add updated recipes for stellar 
evolution and SN explosion to V-Body simulations, even 
though it can be used as a simple and fast stand-alone 
population-synthesis code too. In particular, we merged it 
with the Starlab public software environment (Portegies 
Zwart et al. 2001) and with an upgraded version of HlG- 
PUs code (Capuzzo-Dolcetta et al. (2013); Spera, in prepa¬ 
ration). Thus, the new code can be used for both population 
synthesis studies of compact-object binaries in the field, and 
for investigating the dynamical evolution of compact objects 
in star clusters. The evolution of compact remnants in star 
clusters is of crucial importance, since star clusters are sites 
of intense dynamical processes, which may significantly af¬ 
fect the formation of X-ray binaries (e.g. Blecha et al. 2006; 
Mapelli et al. 2013; Mapelli & Zampieri 2014), as well as 
the formation and merger of double-compact object bina¬ 
ries (e.g. O’Leary et al. 2006; Sadowski et al. 2008; Downing 
et al. 2010, 2011; Ziosi et al. 2014). Furthermore, extreme dy¬ 
namical processes, such as repeated mergers of compact rem¬ 
nants (Miller & Hamilton 2002) and the runaway merger of 
massive objects in star clusters (Portegies Zwart & McMil¬ 
lan 2002), can lead to the formation of intermediate-mass 
BHs (i.e. BHs with mass 10 2 — 10 5 M 0 ). Finally, compact 
remnants are also expected to affect the overall dynamical 
evolution of star clusters (Downing 2012; Sippel et al. 2012; 
Mapelli & Bressan 2013; Trani et al. 2014). 

This paper is organized as follows. In Section 2, we de¬ 
scribe the main features and ingredients of SEVN (including 
stellar evolution and SN models). In Section 3, we discuss 
the outputs of SEVN, with particular attention to the mass 
spectrum and the mass function of NSs and BHs. Further¬ 
more, we compare the results of SEVN with those of other 
population-synthesis codes. In Section 4, we discuss the re¬ 
sults we obtained applying the O’Connor & Ott (2011) and 
Ertl et al. (2015) prescriptions for SN explosion to PAR¬ 
SEC progenitors, at metallicity Z = 0.02. In Section 5, we 
summarize our main results. 


1 SEVN upon request to the authors, through the email 
mario.speraOoapd.inaf.it or mario.speraOlive.it 



The mass spectrum of compact remnants from the PARSEC stellar evolution tracks 3 


2 METHOD 

2.1 Single stellar evolution with PARSEC 

The PARSEC database includes updated and homogeneous 
sets of canonical single stellar evolutionary tracks, from very 
low (M=0.1M@) to very massive (A/=350Mq) stars, and 
from the pre-MS to the beginning of central carbon burn¬ 
ing. The code is thoroughly discussed in Bressan et al. 
(2012, 2013), Chen et al. (2014) and Tang et al. (2014) 
and here we briefly describe its most important character¬ 
istics. The equation of state (EOS) is computed with the 
FreeEOS code 2 (A.W. Irwin). Opacities are computed com¬ 
bining the high-temperature data from the Opacity Project 
At Lawrence Livermore National Laboratory (OPAL) (Igle- 
sias & Rogers 1996) with the low-temperature data from 
the /ESOPUS 3 code (Marigo & Aringer 2009). Conductive 
opacities are included following Itoh et al. (2008). The main 
Hydrogen and Helium burning reactions are included as rec¬ 
ommended in the JINA database (Cyburt et al. 2010) with 
electron screening factors taken from Dewitt et al. (1973) 
and Graboske et al. (1973). Energy losses by electron neu¬ 
trinos are taken from Munakata et al. (1985) and Itoh & 
Kohyama (1983) and Haft et al. (1994). Instability against 
convection is tested by means of the Schwarzschild criterion 
and, where needed, the convective temperature gradient is 
estimated with the mixing-length theory of Bohm-Vitense 
(1958) with a mixing length parameter calibrated on the 
solar model, «mlt = 1.74. The location of the boundary 
of the convective core is estimated in the framework of the 
mixing-length theory, allowing for the penetration of convec¬ 
tive elements into the stable regions Bressan et al. (1981). As 
thoroughly described in Bressan et al. (2013), the main pa¬ 
rameter describing core overshooting is the mean free path 
of convective elements across the border of the unstable re¬ 
gion l c =A c Hp with A c = 0.5, as result of the calibration 
obtained by the analysis of intermediate age clusters (Gi- 
rardi et al. 2009) as well as individual stars (Kamath et al. 
2010; Deheuvels et al. 2010; Torres et al. 2014). Effects of 
stellar rotation have not yet been introduced in PARSEC. 

The reference solar partition of heavy elements is taken 
from Caffau et al. (2011) who revised a few species of the 
Grevesse & Sauval (1998) compilation. According to Caffau 
et al. (2011) compilation, the present-day Sun’s metallicity 
is Z Q = 0.01524. 

While the evolution below M = 12 Mq is computed at 
constant mass, for more massive stars the mass loss rate is 
taken into account combining the mass-loss rates formula¬ 
tions provided by different authors for different evolutionary 
phases, as described in Tang et al. (2014). During the Blue 
Super Giant (BSG) and LBV phases we adopt the maximum 
between the relations provided by Vink et al. (2000, 2001b), 
and that provided by Vink et al. (2011) which includes the 
dependence of the mass-loss rates on the ratio (P) of the 
star luminosity to the corresponding Eddington luminosity. 
In the Red Supergiant (RSG) phases we adopt the mass- 
loss rates by de Jager et al. (1988), R dj, while, in the WR 
phases, we use the Nugis & Lamers (2000) formalism. 

An important effect of the metallicity is its modulation 
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of the mass loss rates. As discussed in Tang et al. (2014) and 
in Chen et al. (2015, in preparation), the dependence of the 
radiation driven mass-loss rates on the metallicity is a strong 
function of T. While, at low values of T, the mass-loss rates 
obey the relation M oc (Z / Zg) 0 ' 85 Mq yr _1 (Vink et al. 
2000, 2001b), with Zg = 0.02 being the average metallic¬ 
ity assumed for Galactic massive stars, at increasing T the 
metallicity dependence becomes weaker, and it disappears 
as r approaches 1 (Grafener & Hamann 2008). Tang et al. 
(2014) show that the metallicity effect can be expressed as 

Mot ( Z/ZaT , ( 1 ) 

with the coefficient a determined from a fit to the pub¬ 
lished relationships by Grafener & Hamann (2008) 

a = 0.85 (T < 2/3) 

a = 2.45 — 2.4 T (2/3 sg T sg 1) (2) 

In the WR phases, PARSEC makes use of the Nugis & 
Lamers (2000) formalism, with its own dependence on the 
stellar metallicity while, during the Red Supergiant (RSG) 
phases the de Jager et al. (1988) rates are re-scaled adopting 
the usual relation M oc {Z/Zg) 0S5 Mq yr _1 . 

With these assumptions for the mass-loss rates, the new 
models of near-solar metallicity can naturally reproduce the 
observed lack of supergiant stars above the Humphreys & 
Davidson (1979) limit. The lack of RSG stars is usually in¬ 
terpreted as a signature of the effects of enhanced mass-loss 
rates when the star enter this region, and this interpretation 
is supported by the presence, around this limit, of LBV stars 
which are known to be characterized by high mass loss rates. 
While, in previous models, the limit was reproduced by 
adopting an “ad-hoc” enhancement of the mass-loss rates, 
in the current models the enhancement is nicely reproduced 
by the boosting of the mass-loss rate when the stars ap¬ 
proach the Eddington limit (Chen et al., in preparation). 
At metallicities lower than solar, the boosting is mitigated 
by the reduction factor introduced by the metallicity depen¬ 
dence. At Z = 0.001, the upper MS widens significantly and 
the more massive stars evolve in the “forbidden” region even 
during the H-burning phase, because of their very large con¬ 
vective cores. They may also ignite and burn central helium 
as “red” super-giant stars. The full set of new evolution¬ 
ary tracks and the corresponding isochrones may be found 
at http://people.sissa.it/~sbressan/parsec.html and 
http: //stev. oapd. inaf . it/cgi-bin/cmd, respectively. 

2.2 SEVN general description 

The coupling between dynamics and stellar evolution, in a 
single code, can be achieved through three alternative ap¬ 
proaches: 

• the first one is based on a “brute force” approach. It 
consists in calling an advanced stellar evolution code (such 
as PARSEC) that calculates the detailed evolution of stellar 
physical parameters step by step, following the time intervals 
imposed by the A-Body dynamics; 

• the second one is based on polynomial fittings that in¬ 
terpolate the fundamental stellar parameters (radius, lumi¬ 
nosity, temperature and chemical composition), as a func¬ 
tion of time, mass and metallicity. Besides being a fast choice 
in terms of computing time, one of the main advantages of 
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using this strategy is that it can be implemented with little 
effort; 

• the third approach consists in using stellar evolution 
isochrones as input files. These isochrones are usually pro¬ 
vided in the form of tables, for a grid of masses and metal- 
licities, and they are read and interpolated by the numerical 
code on the fly. The main advantage of this strategy is that 
it makes the implementation more general. The option to 
change the built-in stellar evolution recipes is left to users, 
who can substitute the input tables, without modifying the 
internal structure of the code or even recompiling it. 

The first approach is highly inefficient because the con¬ 
tinuous calls to advanced stellar evolution codes, inside an 
IV-Body integrator, significantly slows down the overall nu¬ 
merical evolution. To develop SEVN, we chose to follow 
the second aforementioned approach (usage of stellar evo¬ 
lution isochrones in tabular form). SEVN can work as a 
stand-alone code (for fast population synthesis studies in 
the field), and can be linked to a large variety of IV-Body 
codes, without suffering a performance penalty. In particu¬ 
lar, we merged SEVN with an updated version of the di¬ 
rect IV-Body code HlGPUs 4 (Capuzzo-Dolcetta et al. 2013; 
Spera, in preparation) as well as in the Starlab software 
environment 5 (Portegies Zwart et al. 2001), and it can also 
be included in the Astrophysical Multipurpose Software En¬ 
vironment (AMUSE 6 , Pelupessy et al. 2013). 

In this paper, we focus our attention on our implemen¬ 
tation of SEVN in Starlab since it already includes both 
an V-Body integrator (called Kira) and a binary evolution 
module (SeBa). In particular, we updated a version of SeBa 
that had been previously modified by Mapelli et al. (2013), 
who included metallicity dependent stellar winds (Hurley 
et al. 2000) and prescriptions for the mass loss by MS stars 
(Vink et al. 2001a). While we left the dynamical integra¬ 
tion part untouched, we rearranged SeBa by adding stellar 
isochrone tables, at different metallicity, and by forcing the 
software to use them as input files. In this way, we have 
hidden the default implementation without making radi¬ 
cal changes to the code structure. In the current version 
of SEVN, we use the PARSEC data to get the physical 
parameters of the stars for all evolutionary stages but the 
thermally-pulsating AGB phase (TP-AGB). In fact, the evo¬ 
lution and lifetimes of TP-AGB stars suffer from significant 
uncertainties and a thorough calibration of the latter phase 
is still underway (Marigo et al. 2013; Rosenfield et al. 2014). 
At present, we use the built-in SeBa super_giant class to 
follow the evolution of the stars in this stage. Moreover, 
according to the PARSEC recipes, all stars with an ini¬ 
tial mass Mzams < M up (with M up = 7M 0 ) undergo the 
AGB phase. In particular, at the end of their lives, stars 
of mass Mzams ^ M up will explode as SNe leaving NSs or 
BHs as compact remnants, while stars with Mzams < M up 
will evolve through the AGB phase, quickly losing their en¬ 
velopes, until a WD is formed. More technical details about 
the SEVN implementation can be found in Appendix A. 


4 http://astrowww.phys.uniromal.it/dolcetta/HPCcodes/ 
HiGPUs.html 

1 http://www.sns.ias.edu/~starlab/ 

6 http://amusecode.org/wiki 


2.3 Prescriptions for the formation of compact 
remnants 

The default recipes implemented in the SeBa module pre¬ 
dict the formation of a white dwarf (WD) if the final core 
mass is less than the Chandrasekhar mass (1.4 M 0 ), a NS 
or a BH if the core mass is greater than 1.4 M 0 . In our im¬ 
plementation of SEVN in SeBa, we leave the recipes for 
the formation of WDs unchanged, but we change the way to 
form NSs and BHs. 

The default version of SeBa distinguishes between NSs 
and BHs by inspecting the final mass of the core: if it is 
larger than the Chandrasekhar mass (1.4 M 0 ) and, at the 
same time, the initial mass of the star is AFzams < 25 M 0 , 
a NS is formed. If Mzams ^ 25 M 0 or if the final carbon- 
oxygen (CO) core mass (Mco) is such that Mco ^ 5M 0 , 
the star ends its life forming a BH 7 . To determine the BH 
mass, SeBa assumes that, initially, a fixed amount of the 
CO core mass collapses, forming a proto-compact object of 
mass Mproto = 3M 0 (Fryer & Kalogera 2001). The amount 
of fallback material, A/fb, is determined by comparing the 
binding energies of the hydrogen (H), helium (He) and CO 
shells with the SN explosion energy. The final mass of the 
compact object is given by Mbh = M prot o + Mfb. 

In SEVN, we substituted the default SEBA treatment 
of SNe with the following new recipes. We implemented the 
three models described in details by Fryer et al. (2012): (i) 
the model implemented in the StarTrack population syn¬ 
thesis code (see Belczynski et al. 2008, 2010), (ii) the rapid 
supernova model, and (iii) the delayed supernova model. The 
main difference between the last two explosion mechanisms 
is the time-scale over which the explosion occurs: < 250 
ms after the bounce for the rapid model, > 0.5 s for the 
delayed mechanism (for the details see, for example, Bethe 
1990). A common feature of these models is that they de¬ 
pend only on the final characteristics of the star, by means 
of the final CO core mass (Mco) and of the final mass of 
the star (Mfi n ). Appendix B summarizes the main features 
of the Fryer et al. (2012) SN explosion recipes. We recall 
that the Fryer et al. (2012) methods are general prescrip¬ 
tions for the formation of compact remnants, and do not 
distinguish, a priori, between NSs and BHs. In SEVN, we 
assume that all the remnants with masses Af rem < 3.0 M 0 
are NSs, and that the objects with masses M rem is 3.0 M 0 
are BHs, according to the maximum mass of a NS indicated 
by the Tolman-Oppenheimer-Volkoff limit (Oppenheimer & 
Volkoff 1939). While the Fryer et al. (2012) models are ex¬ 
tremely simple to implement in a population-synthesis code, 
it has been recently suggested that the dependence of the 
mass of the compact remnant on A/fi n or Mco might be 
significantly more complex (Ugliano et al. 2012; O’Connor 

6 Ott 2011; Suklibold & Woosley 2014; Janka 2012; Smartt 
2015). The internal structure of stars, at core-collapse stage, 
may exhibit significant differences, leading to deep changes 
on the physical parameters of compact remnants, even if the 
progenitors are very close in terms of AFzams or Mco- As 
a consequence, a one to one relation between the mass of 

7 In SeBa, the limits 25 M 0 and 5 M 0 are the default val¬ 
ues of two parameters called super _giant2blackJiole and 
C0core2black_hole, respectively. The user can adjust them at 
choice. 
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the compact remnant and, e.g., Mco could be inadequate 
to discriminate between SNe (formation of a NS) and failed 
SNe (direct collapse to a BH). The critical parameter to 
distinguish between SNe and failed SNe might the compact¬ 
ness of stellar cores at the pre-SN stage (O’Connor & Ott 
2011; Ugliano et al. 2012; Sukhbold & Woosley 2014). Al¬ 
ternatively one may use an equivalent criterion based on 
the two-parameters M 4 , representing the enclosed mass at 
a dimensionless entropy per nucleon s = 4, and /m, that is 
the mass gradient at the same location (Ertl et al. 2015). 
In order to fulfil these recent advances of the SN explosion 
models and to test their impact on the mass spectrum of 
compact remnants, we have implemented in SEVN these 
two additional SN explosion recipes, namely the criterion 
based on the compactness of stellar cores (O’Connor & Ott 
2011; Ugliano et al. 2012; Sukhbold & Woosley 2014) and the 
criterion based on M 4 and /m (Ertl et al. 2015). We present 
here only the results for Z = 0.02, because the results at 
metallicities lower than solar are still under investigation. 

I 11 SEVN, we set the delayed SN model as default SN 
explosion mechanism, but the user can choose one of the 
aforementioned mechanisms by modifying the input param¬ 
eter file. Only for the SEVN implementation in Starlab, 
we also leave the choice to use the SeBa built-in recipes 8 . 

Furthermore, the aforementioned models do not ac¬ 
count for the possibility that the progenitor undergoes a 
pair-instability SN (e.g. Woosley et al. 2002). In SEVN, 
we add the option to activate pair-instability SNe, when 
the Helium core mass (after the He core burning phase) is 
60 ^ Mne/ M@ ^ 133. For this range of He core masses, the 
star does not leave any remnant, while it directly collapses 
to BH for larger masses. In the following Section, we show 
models that do not undergo pair-instability SNe. 

When a compact remnant is formed, it also receives a 
velocity kick, Wkick, due to the asymmetries that can occur 
during the collapse process. In SEVN, we determine the ab¬ 
solute value of the kick using the three dimensional velocity 
distribution of the pulsars observed in our galaxy. For de¬ 
tails, we refer to Hobbs et al. (2005), who studied the proper 
motions of 233 pulsars, obtaining a Maxwellian fit for their 
velocity distribution, with a one dimensional variance equal 
to 256 km/s. The direction of the kick is randomly cho¬ 
sen. Furthermore, following the prescriptions given in Fryer 
et al. ( 2012 ), we also included the dependence of the ve¬ 
locity kick on the amount of mass that falls back onto the 
proto-compact object. Specifically, the actual value of the 
kick imparted to a compact remnant, I 4 ick, is given by 

Vkick = (1 — /fb) Wkick- (3) 

Thus, a BH that forms via direct collapse (/fb = 1) does 
not receive a velocity kick, while full kicks are assigned to 
compact remnants formed with no fallback. Another possi¬ 
ble treatment for BH kick velocities, is to assume that BHs 
follow the same distribution of Wkick as NSs, but normalized 
to (Mns)/Mbh (where (Mns) is the average NS mass), to 

8 The first line of the file input_param.txt determines the SN 
explosion model that will be adopted throughout the numerical 
simulation. It can be delayed, startrack, rapid, compactness 
or twoparameters. In the implementation of SEVN in Starlab, 
this line can be even default if we want to use the SeBa built-in 
recipes for SN explosion. 


Z = 2.0E-2 A Final mass PARSEC □ Final mass SSE 



Figure 1. Temporal evolution of the stellar mass for different 
Mzams an( l for metallicity Z = 2 .Ox 10 2 , according to PARSEC 
and SSE. Solid (dashed) line: evolution of a star with Mzams = 
10 Mq with PARSEC (SSE); dotted (dash-dotted) line: evolution 
of a star with Mzams — 30 Mq with PARSEC (SSE); dash- 
double dotted (short dashed) line: evolution of the mass of a star 
with Mzams = 60 Mq with PARSEC (SSE); short dotted (short 
dasli-dotted) line: evolution of the mass of a star with Mzams = 
110 Mq with PARSEC (SSE). Open triangles and open squares 
mark the final point of each curve obtained using PARSEC and 
SSE, respectively. 

ensure momentum conservation. We leave this second option 
in Starlab, even if we set the former treatment as default. 


3 RESULTS 

In this section, we discuss the effects of metallicity on the 
stellar mass loss rate, on the CO core, on the formation 
of compact remnants, and on the mass function of NSs and 
BHs, as we found using our new tool SEVN. We also discuss 
the main differences between SEVN and other population 
synthesis codes, in terms of mass spectrum of compact rem¬ 
nants. In particular, we compare the results of SEVN with 
those of SSE (Hurley et al. 2000), of Starlab (Portegies 
Zwart et al. 2001), and of the version of Starlab modified 
by Mapelli et al. (2013) (hereafter referred as StarlabMM). 
In particular, SSE is a stellar evolution tool that has al¬ 
ready been linked to the NBODYx family of WBody codes 
(see e.g. Aarseth (1999) and Nitadori & Aarseth (2012)) and 
it also implements recipes for metallicity-dependent stellar 
winds. Moreover, SSE adopts the SN explosion recipes de¬ 
scribed in Belczynski et al. (2002). 

3.1 Mass loss by stellar winds 

Figures 1, 2 and 3 show the temporal evolution of stellar 
mass at Z = 2 x 10 -2 , 2 x 10 ~ 8 and 2 x 10 -4 , respectively, 
for four selected ZAMS masses between 10 and 110 Mq. 
The evolution of the stellar mass predicted by PARSEC is 
compared with that implemented in SSE. At lower ZAMS 
masses (Mzams ;$ 10 Mq) the behaviour of PARSEC and 
SSE is almost indistinguishable. 
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Z = 2.0E-3 A Final mass PARSEC □ Final mass SSE 



Figure 2. Same as Fig. 1, but for Z = 2.0 X 10 3 . 


A Final mass parsec O Final mass SSE 



Time (Myr) 


Figure 3. Same as Fig. 1, but for Z = 2.0 X 10 4 . 

For larger masses, there is no significant difference for 
most of the star’s life, but there is a significant difference 
in the final masses M& n , especially at low metallicity. The 
differences in Ma n are about 80% of Mzams for stars with 
Mzams > 60 M© at Z = 2 x 10 -4 . The reason of these 
differences is the treatment of stellar winds, especially in the 
late-MS, LBV and WR stages (see Section 2.1 and Bressan 
et al. 2012; Tang et al. 2014 for details). 

3.2 Final mass (Man) and CO core mass (Meo) 

The SN explosion mechanisms discussed by Fryer et al. 
(2012), and implemented in SEVN, depend on the final mass 
of the star, M& n , and on its final CO mass, Mco ( see equa¬ 
tions B2, B5 and B8). Since both Mfi n and Mco depend on 
the initial mass of the star, Mzams, and on its metallicity, 
Z, thus also M r em will depend on Mzams and Z. This im¬ 
plies that the mass spectrum of compact remnants strongly 
depends on the prescriptions adopted to evolve the star until 
its pre-SN stage. 



Figure 4. Final mass of the stars as a function of their initial 

mass, for different values of metallicity 1.0 X 10 -4 . Z 2.0 X 

10 -2 . Top solid line: Z = 1.0 X 10 -4 ; dashed line: Z = 2.0 X 10 -4 ; 
dotted line: Z = 5.0 X 10 4 ; dash-dotted line: Z = 1.0 X 10 -3 ; 
dash-double dotted line: Z = 2.0 X 10 -3 ; short dashed line: Z = 
4.0 X 10 -3 ; short dotted line: Z = 6.0 X 10 -3 ; short dash-dotted 
line: 1.0 X 10 -2 ; bottom solid line: Z = 2.0 X 10 -2 . 



Figure 5. Final mass of the CO core as a function of the initial 
mass of the star, for different values of metallicity 1.0 X 10 -4 T 
Z . 2.0 X 10 -2 . Line types are the same as in Fig. 4. 

Figure 4 shows the trend of Mfi n as a function of 
Mzams, for different values of the metallicity. Figure 4 re¬ 
flects the fact that metal-poor stars are subject to weaker 
stellar winds throughout their evolution. In fact, Ata n is al¬ 
ways smaller than ~ 25 M@ at Z = 2.0 X 10~ 2 , while Mfi n « 
Mzams at Z < 2.0 x 10 -4 . The curves for Z < 2.0 x 10 -4 
are well approximated by a simple linear relation 

Mfin (Mzams) = 0.9519Mzams + 1-45. (4) 

In Fig. 5, we show Mco as a function of Mzams, 
for different values of metallicity. As expected, the final 
CO mass scales inversely with metallicity: the maximum 
value of Mco ranges between ~ 20 Mq and ~ 65 Mg, for 
1.0 x 10 -4 ^ Z ^ 2.0 x 10 -2 . It is interesting to note that, 
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Figure 6. Mass of the final compact remnant (M rem ) as a func¬ 
tion of the initial mass of the star, for various metallicities. The 
curves have been obtained using SEVN and the delayed SN 
model. Line types are the same as in Fig. 4. 



Figure 7. Mass of the compact remnant as a function of the 
final CO core mass of the progenitor, for different metallicities. 
The curves have been obtained using SEVN and the delayed SN 
model. Line types are the same as in Fig. 4. 


for Z ^ 1.0 x 10 3 , the curves of Fig. 5 become approxi¬ 
mately independent of Z, and can be expressed as 


'0.3403A7 Z ams - 2.064 


Mqo (Mzams) = 


if Mzams < 27 M@ 
0.4670M ZA ms - 5.47 


if Mzams ^ 27 M 0 


(5) 


At present, since PARSEC does not include the TP-AGB 
stellar evolution phase, equation 5 holds for Mzams > 
M up = 7M@ (see Sec. 2.2). 


3.3 The mass spectrum of compact remnants 

In Fig. 6, we show the mass spectrum of compact rem¬ 
nants as a function of the ZAMS mass of their progenitors, 
for different values of metallicity. To obtain the curves in 
Fig. 6, we used the delayed supernova model, chosen as 
the default explosion mechanism in SEVN. As expected, 
in Fig. 6, we notice that the lower the metallicity is, the 
higher the mass of the heaviest compact remnant; in par¬ 
ticular, M rem ranges from ~ 25 Mq at Z = 2.0 x 10~ 2 to 
~ 135 M 0 at Z = 1.0 x 10 -4 . For Z < 2.0 x 10~ 4 and 
7 Mq = M up ^ Mzams $ 150 M 0 , simple fitting formulas 
can be derived for M rem (Mzams), by substituting the best 
fit curves for M Rn (A/ Z ams, Z) and Meo (A7 Z ams, Z) (Eqs. 4 
and 5, respectively) in the formulas of the delayed explosion 


mechanism (Eq. B8): 

"l.4M 0 

if M up < Mzams 7; 13 M 0 


0.170Mzams — 0.882 

if 13 M 0 < Mzams 1$ 16 M 0 


M r , 


(0.041M| ams — 0.673M|ams+ 

+2.18A7 Z ams + 0.361)/ (0.952M Z ams + 0.15) 
if 16 M 0 < Mzams 1$ 27 M 0 


(0.0563M| ams — 1.10M| ams + 

V2.49A7zams + 0.318)/ (0.952M Z ams + 0.15) 
if 27 M 0 < Mzams 1$ 36 M 0 


0.952M Z ams + 1-45 


if A7 Z ams ^ 36 M 0 . 

( 6 ) 

A general fitting formula for M rem , as a function of Mzams 
and Z, and that holds for every metallicity, is provided in 
Appendix C. 

Figure 7 shows the value of M rem as a function of Mco, 
for different metallicities. It is worth noting that, for ev¬ 
ery metallicity, A/ rem lies approximately between M rem ,u P = 
1.85Mco + 11.9 and A/ rem ,down = 1.22Mco + 1-06 (see Ap¬ 
pendix C for the details). 


3.4 Comparison of different supernova explosion 
models 

In Fig. 8, we show the mass of the remnants as a func¬ 
tion of Mzams, for different SN recipes, at fixed metallicity 
Z = 2.0 x 10 -2 , in order to compare the various SN mod¬ 
els implemented in SEVN. In this figure, we also show the 
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Figure 8. Mass of compact remnants as a function of MzaMS at 
Z = 2.0 X 10 -2 , derived with SEVN using different models of SN 
explosion. Dash-double dotted line: StarTrack SN recipes; solid 
line: delayed SN model; dashed line: rapid SN model; dotted line: 
Starlab prescriptions. 


results obtained using the SeBa (Starlab) built-in mod- 
els(see Portegies Zwart et al. 2001 for details). 

Figure 8 shows that all recipes produce approximately 
the same remnant mass spectrum, for Mzams > 50 Mg. The 
bottom panel of Fig. 9 is a zoom of Fig. 8 in the region of 
25 Mq ^ Mzams $ 50 Mq. From Fig. 9 we notice that the 
StarTrack SN recipes produce, on average, more massive 
BHs (with mass between ~ 12 Mq and ~ 18 Mq) in the in¬ 
terval 28 Mq < Mzams ;$ 50 M 0 . This is due to the fact 
that StarTrack predicts the formation of compact rem¬ 
nants via direct collapse if Mco ^ 7.6 M 0 (see equation 
B2), condition that occurs for A/zams > 28M 0 (see Fig. 
5). The other models do not predict direct collapse in this 
interval of A/zams, and produce lighter BHs with masses 
between ~ 6M 0 and ~ 13 M 0 . 

The abrupt step of the rapid SN model, for 24 M 0 < 
Mzams ;$ 26 M 0 , corresponds to the process of direct col¬ 
lapse that takes place for 6 M 0 ^ Mco ^ 7 M 0 in this model 
(see equation B5). 

In the range 14 M 0 ^ Mzams ^ 24 M 0 (see upper panel 
of Fig. 9), the delayed model predicts a higher amount of 
fallback than the other models. In fact, the delayed mecha¬ 
nism forms compact objects with masses between ~ 2.0 M 0 
and ~ 6.0 M 0 , while the other models form remnants with 
masses only up to ~ 2 M 0 (see upper panel of Fig. 9). Using 
the StarTrack prescriptions, it is possible to form rem¬ 
nants with masses > 3.0 M 0 , but only for Mzams •Sj 22 M 0 . 
Finally, using the SN model implemented in SeBa and the 
rapid SN model, we find a paucity of remnants with masses 
between ~ 2 M 0 and ~ 6 M 0 with the result of having a 
marked gap between the heaviest NS and the lightest BH. 

In Figs. 10 and 11, we show the mass spectrum of 
BHs and NSs obtained for different explosion models for 
Z = 2.0 x 10~ 3 and Z = 2.0 x 10~ 4 , respectively. At 
Z = 2.0 x 10~ 3 (Z = 2.0 x 10~ 4 ), the maximum BH mass is 
~ 60 M 0 (~ 130 M 0 ), regardless of the SN explosion mech¬ 
anism. The main remarkable features of Figures 10 and 11 
are the following: 




Figure 9. Two details of Fig. 8: the top panel shows the range 
8M 0 < Mzams $ 25 M 0 , while the bottom panel refers to the 
interval 25 Mq < Mzams $ 50 M 0 . In the top panel, the hor¬ 
izontal dashed line marks the transition between NSs and BHs. 
The thick, semi-transparent line in the bottom panel highlights 
the intervals in which direct collapse occurs. The other lines are 
the same as in Fig. 8. 


(i) the StarTrack models produce heavier compact 
remnants for 25 M 0 ^ Mzams ^ 35 M 0 ; 

(ii) the rapid SN model exhibits an abrupt step for 
24 Mq ^ A/zams ^ 26 M 0 ; 

(iii) for Mzams > 35 M 0 , the mass spectra obtained with 
the models of Fryer et al. (2012) become indistinguishable 
(all of them predict direct collapse); 

(iv) except for the delayed model, we obtain a paucity of 
remnants with masses between ~ 2M 0 and ~ 6M 0 ; 

(v) The SeBa built-in SN explosion model predicts direct 
collapse for Mzams > 45 M 0 (Mzams > 40 M 0 ) at Z = 
2.0 x 10 -3 (Z = 2.0 x 10 -4 ). 

In Section 4, we extend this comparison to more sophis¬ 
ticated models of SN explosion (based on the compactness 
of the stellar core at pre-SN stage, and on the dimensionless 
entropy per nucleon at pre-SN stage). 
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Figure 10. Mass of compact remnants as a function of A'// a vis 
at Z = 2.0 X 10~ 3 , derived with SEVN using different models of 
SN explosion. Line types are the same as in Fig. 9. 



Figure 11. The same as Fig. 10, but for Z = 2.0 X 10 4 . 

3.5 Comparisons with other stellar evolution tools 

Figure 12 shows the mass spectrum of compact remnants, 
at Z = 2.0 x 10 -2 , obtained using SEVN, with the de¬ 
layed supernova explosion model, in comparison with the 
results of Starlab v4.4.4 (default SeBa stellar evolution 
module, Portegies Zwart et al. 2001, hereafter simply Star- 
lab), StarlabMM (Mapelli et al. 2013) and SSE. The max¬ 
imum BH mass we obtain, at Z = 2. Ox 10 -2 , using SEVN, is 
~ 25 Mg, while using Starlab this value is slightly higher 
(~ 28 Mg). In SEVN, the stars with Mzams — 100 Mg 
form the heaviest BHs, while, using Starlab, the most mas¬ 
sive remnants derive from stars with 85 Mg ^ Mzams 5) 
150 Mg . StarlabMM produces BHs with masses up to 
~ 23 Mg. It is interesting to point out that the recipes imple¬ 
mented in Starlab produce a paucity of compact remnants 
with masses between ~ 2 Mg and ~ 5 Mg . This gap derives 
from the assumption that BHs form only if Meo ^ 5 Mg, 
otherwise, NSs with masses between ~ 1.2 Mg and ~ 1.6 Mg 
are formed. If we use the SSE package, the maximum mass 



Figure 12. Mass of compact remnants as a function of Mzams a t 
Z = 2.0 X 10 -2 , derived with different codes: SEVN (solid line), 
SSE (dotted line), StarlabMM (dashed line), Starlab v 4.4.4 
(dash-dotted line). The semi-transparent line highlights the in¬ 
tervals in which direct collapse takes place. For SEVN, we used 
the delayed SN mechanism. StarlabMM is the modified version 
of Starlab described in Mapelli et al. (2013), while Starlab 
v4.4.4 is the standard version of Starlab (ver. 4.4.4). 

of compact remnants is ~ 13Mg. It is also important to 
stress that, for 17Mg < Mzams 1$ 40 Mg, the delayed ex¬ 
plosion model implemented in SEVN creates more massive 
compact remnants than the other models. 

Figures 13 and 14 show the mass spectrum of com¬ 
pact remnants at Z = 2.0 X 10~ 3 and Z = 2.0 x 10 —4 , 
respectively. The results of Starlab are not shown in Fig¬ 
ures 13 and 14, because Starlab does not include metal- 
licity dependent stellar winds. At Z = 2.0 x 10~ 3 and for 
Mzams ^ 30 Mg, SEVN (with the PARSEC evolution¬ 
ary tables) produces significantly heavier BHs than Star¬ 
labMM and SSE (see Fig. 13). In particular, the maxi¬ 
mum BH mass obtained using SEVN is ~ 60 Mg, while 
this value is ~ 40 Mg and ~ 20 Mg in the case of Star¬ 
labMM and SSE, respectively. We also stress that, while for 
SEVN and StarlabMM the heaviest BH comes from the 
death of the most massive star (that is Mzams = 150 Mg), 
in the case of SSE BHs of ~ 20 Mg form from stars with 
25 Mg < Mzams ;$ 30 Mg only. The abrupt step observed 
for Mzams — 100 Mg, in the StarlabMM curve represents 
the transition between partial fallback and direct collapse 
(occurring at Mfi n ^ 40 Mg, see Mapelli et al. 2013 for de¬ 
tails), while that at Mzams — 25 Mg reflects the transition 
from NSs to BHs. It is worth noting that, for Z = 2.0 x 10 -3 
and Mzams 30 Mg, the SSE model produces more mas¬ 
sive compact remnants than the other models. 

Similar considerations hold for Fig. 14, which is the 
same as Figs. 12 and 13 but at Z = 2.0 x 10~ 4 . SEVN 
predicts BHs masses up to ~ 120 Mg, StarlabMM cre¬ 
ates BHs of maximum mass ~ 80 Mg, while the SSE pre¬ 
scriptions do not go beyond ~ 25 Mg. Also in this case, as 
observed at Z = 2.0 x 10 —3 , the SSE recipes predict the 
formation of more massive compact remnants in the range 
20 Mg < AfzAMS )$ 30 Mg. 

Tables 1, 2 and 3 report the values of Mzams, Mfi n and 
Mco, corresponding to the transition between the forma- 
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Figure 13. Mass of compact remnants as a function of Mzams 
at Z = 2.0 X 10 —3 , derived with different codes: SEVN (solid 
line), SSE (dotted line), StarlabMM (dashed line). The semi¬ 
transparent line highlights the intervals in which direct collapse 
takes place. For SEVN, we used the delayed SN mechanism. 



Figure 14. The same as Fig. 13 but for Z = 2.0 X 10 4 . 


Table 1 . Values of Mzams > Mg n and Mco that correspond 
to the transition between the formation of a NS and that of a 
BH, and maximum BH mass (M™^ x ), for three different codes: 
SEVN, StarlabMM and SSE. D: delayed model; R: rapid model; 
S: StarTrack prescriptions. Results for Z = 2.0 X 10“ 2 . 



D 

SEVN 

R 

S 

StarlabMM 

SSE 

Mzams 

18.8 

23.9 

21.8 

23.0 

20.7 

Mfl n 

16.0 

19.6 

18.7 

17.3 

7.3 

Mco 

4.1 

6.0 

5.1 

5.0 

5.3 

m bh x 

25.0 

25.0 

25.0 

23.0 

12.0 


Table 2. 

Same as 

Tab.l but for Z = 

2.0 x 10“ 3 . 



D 

SEVN 

R 

S 

StarlabMM 

SSE 

Mzams 

18.2 

23.6 

21.3 

23.0 

18.8 

Mfi n 

17.9 

23.2 

21.1 

17.6 

16.4 

Mco 

4.1 

6.0 

5.1 

5.0 

5.1 

m bh x 

58.0 

58.0 

58.0 

39.0 

19.0 


Table 3. Same as 

Tab.l but for Z = 

2.0 x 10“ 4 . 



D 

SEVN 

R 

S 

StarlabMM 

SSE 

Mzams 

18.2 

23.1 

21.3 

23.0 

18.0 

Mfln 

18.0 

23.1 

21.3 

17.7 

17.0 

Mqo 

4.1 

6.0 

5.1 

5.0 

5.1 

Mbh x 

130.0 

130.0 

130.0 

83.0 

26.0 


tion of a NS and a BH, at Z = 2 X 10 2 , 2 x 10 -3 and 
2 x 10 —4 , respectively. The results obtained using SEVN, 
StarlabMM and SSE are compared in the Tables. We no¬ 
tice that the transition value of Mco does not depend on 
metallicity and it ranges from ~ 4.0 M 0 (delayed model of 
SEVN) to ~ 6.0M 0 (rapid model of SEVN). The transi¬ 
tion values of Mzams and Ma n show a weak dependence on 
metallicity for a given code. Mzams goes form ~ 18 Mq (de¬ 
layed model of SEVN at low metallicity) to ~ 24 M 0 (rapid 
model of SEVN at Z = 2 x 10 -2 ), while Mfi n ranges form 
~ 7M 0 (SSE at Z = 2 x 10 -2 ) to 23 M 0 (rapid model of 
SEVN at low metallicity). In the last row of tables 1, 2 and 
3, we also report the maximum compact remnant mass. As 
we have already shown in this section, for the maximum BH 
mass we get huge differences between the considered codes. 
This is due to the different stellar evolution recipes adopted 
in PARSEC, SSE and StarlabMM, especially for metal- 
poor stars. 

3.6 The mass distribution of compact remnants 

In this section, we derive the mass function of compact 
remnants (NSs and BHs) that form in a stellar population 
following the Kroupa initial mass function (IMF, Kroupa 
2001). The Kroupa IMF scales as AN/Am oc m~ a , with 
a = 1.3 (2.3) for m < 0.5 M 0 (> 0.5 M 0 ). We assume 
a minimum mass m m i n = 0.1 M 0 and a maximum mass 
m m ax = 150 M 0 . We consider three different metallicities 
(Z = 2 x 10 —2 , 2 x 10 -3 and 2 x 10~ 4 ). For each metallic¬ 
ity, we generate 2.5 x 10 6 MS stars, with mass distributed 
according to the Kroupa IMF, and we evolve them with 
SEVN. For each case, we do three realizations: one with the 
delayed SN model, one with the rapid SN model, and one 
with the Startrack recipes for compact remnants. Moreover, 
we also compare SEVN (with the delayed SN recipe) with 
StarlabMM and with SSE. Table 4 lists the properties of 
the different realizations. 

Figure 15 shows the mass distribution of compact rem¬ 
nants obtained for runs Z1D, Z1R and Z1S (see Table 4). 
These stellar populations have Z = 2.0 x 10 -2 and are 
evolved using SEVN, with the PARSEC stellar evolution 
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Table 4. General properties of the stellar populations used as 
test-case to study the mass distribution of NSs and BHs. 


Run 

Z 


SN recipe 

Code 

Z1D 

2 x 

10~ 2 

delayed SN 

SEVN 

Z2D 

2 x 

10“ 3 

delayed SN 

SEVN 

Z3D 

2 x 

10- 4 

delayed SN 

SEVN 

Z1R 

2 x 

10" 2 

rapid SN 

SEVN 

Z2R 

2 x 

10 -3 

rapid SN 

SEVN 

Z3R 

2 x 

10 -4 

rapid SN 

SEVN 

Z1S 

2 x 

10 -2 

Startrack SN 

SEVN 

Z2S 

2 x 

10" 3 

Startrack SN 

SEVN 

Z3S 

2 x 

10" 4 

Startrack SN 

SEVN 

Z1MM 

2 x 

10 -2 

Mapelli et al. (2013) 

StarlabMM 

Z2MM 

2 x 

10“ 3 

Mapelli et al. (2013) 

StarlabMM 

Z3MM 

2 x 

10" 4 

Mapelli et al. (2013) 

StarlabMM 

Z1SSE 

2 x 

10" 2 

Hurley et al. (2000) 

SSE 

Z2SSE 

2 x 

10“ 3 

Hurley et al. (2000) 

SSE 

Z3SSE 

2 x 

10 -4 

Hurley et al. (2000) 

SSE 


We generated and evolved 2.5 X 10® stars in each of these runs. 


Total number of stars n = 2.5xio 6 



Figure 15. Fraction of compact remnants, normalized to the 
total number of stars (JV = 2.5 X 10®) that initially follow a 
Kroupa IMF. Solid line with open triangles: SEVN with delayed 
SN model (Z1D); dash-double dotted line with circles: SEVN 
with rapid SN model (Z1R); dash-dotted line with open circles: 
SEVN with StarTrack recipes (Z1S). The vertical dashed line 
at Mrem = 3 Mg distinguishes NSs from BHs. The curves have 
been obtained for Z = 2.0 X 10 -2 . 

prescriptions and with different SN models (delayed SN 
model, rapid SN model and Startrack recipes for run Z1D, 
Z1R and Z1S, respectively). Both the delayed and rapid 
models predict a peak of BHs with mass ~ 10 Mg at 
Z = 2.0 X 10 -2 , while this peak is shifted to ~ 13Mg in 
the StarTrack prescriptions. The reason for these peaks 
can be understood from Fig. 9: for example, in the delayed 
model, BHs of mass 9 Mg < Mbh < 11 Mg can form from 
a wide range of stars (those with 26 Mg < Mzams 1$ 28 Mg 
and with 35 Mg < Mzams iS 44 Mg). 

Fig. 15 also shows that the rapid SN model predicts 
almost no remnants with mass between ~ 2 Mg and ~ 5 Mg. 


Total number of stars n = 2.5xio 6 



Figure 16. Same as Fig. 15 but for Z = 2.0 X 10 3 . 


Total number of stars n = 2 . 5 x 10 s 



M (Ml 

rem ' O' 


Figure 17. Same as Fig. 15 but for Z = 2.0 X 10 4 . 

This agrees with current observations, which suggest a gap 
between the maximum NS mass and the minimum BH mass 
(Ozel et al. 2010) 9 . Figures 16 and 17 are the same as Fig. 15, 
but for Z = 2.0 x 10~ 3 (runs Z2D, Z2R, Z2S) and Z = 2.0 x 
10 -4 (runs Z3D, Z3R, Z3S), respectively. In these Figures, 
the peak of BH mass distribution is at ~ 35 — 40 Mg. 

The mass distribution for NSs peaks at 1.3 — 1.6 Mg for 
all the SEVN models, almost independently of metallicity. 
A relevant difference between the models is that the delayed 
SN model forms a not negligible number of NSs with masses 
between 2 Mg and 3 Mg while, for the other SN explosion 

9 Whether the presence of this gap is physical or simply due to 
selection biases is still unclear (Farr et al. 2011; Kochanek 2014; 
Ugliano et al. 2012), especially after the recent estimation of the 
BH mass of the X-ray source SWIFT J1753.5-0127, which seems 
to fall right into this gap (Neustroev et al. 2014). 
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Total number of stars n = 2.5 xio 6 



Figure 18. Mass function of NSs and BHs, normalized to the 
total number of stars (N = 2.5 X 10®) that initially follow a 
Kroupa IMF, obtained using three different codes: SEVN, Star- 
labMM and SSE. Solid line with open triangles: SEVN with the 
delayed SN model (Z1D); dash-double dotted line with circles: 
StarlabMM (Z1MM); dash-dotted line with open circles: SSE 
(Z1SSE). The vertical dashed line at M rem = 3 Mq, distinguishes 
NSs from BHs. The curves have been obtained for Z = 2.0 X 10~ 2 . 



Figure 19. Same as Fig. 18 but for Z = 2.0 X 10 3 . 


mechanisms, the vast majority of NSs have a mass below 
2 Mq. 

Figure 18 compares the mass distribution of compact 
remnants in runs Z1D, Z1MM and Z1SSE (i.e. the same 
stellar population at metallicity Z = 2.0 X 10~ 2 , run with 
SEVN, StarlabMM and SSE, respectively). Run Z1MM 
(StarlabMM) agrees with run Z1R (SEVN with the rapid 
SN mechanism) to reproduce the dearth of compact rem¬ 
nants with mass between ~ 2Mq and ~ 5Mq. The major¬ 
ity of BHs formed in run Z1SSE (SSE) have mass between 
~ 5Mq and ~ 12 Mq. SSE does not produce compact rem- 


Total number of stars n = 2.5xio 6 



M ( M ) 

rem ' Q ' 

Figure 20. Same as Fig. 18 but for Z = 2.0 X 10 -4 . 


Table 5. Fraction of BHs, normalized to total number of stars, 
obtained from SEVN, adopting different SN explosion models, 
and from StarlabMM and SSE. D: delayed model; R: rapid 
model; S: StarTrack prescriptions. 


z 


SEVN 


StarlabMM 

SSE 


D 

R 

S 



2.0 x 10“ 4 

2.38 

1.72 

1.94 

1.72 

2.40 

2.0 x 10“ 3 

2.40 

1.66 

1.92 

1.72 

2.28 

2.0 x 10“ 2 

2.26 

1.62 

1.86 

1.72 

2.02 


The values are normalized to 10 3 . 


nants with mass > 13 Mq. As to NSs, SSE produces more 
NSs with masses between ~ 1.5 Mq and 2.5 Mq than the 
other codes, while StarlabMM does not form NSs with 
mass > 1.5 Mq. 

Figures 19 and 20 are the same as Fig. 18 but for 
Z = 2.0 x 10~ 3 and Z = 2.0 x 10 -4 , respectively. At low 
metallicities, SEVN produces heavier BHs than both Star¬ 
labMM and SSE. The majority of BHs in both run Z2SSE 
and Z3SSE have mass ~ 10 — 20 Mq, while the BH mass 
in both run Z2D and Z3D peaks at about ~ 40 Mq. In 
run Z2D (run Z3D) the distribution of BH masses extends 
up to ~ 60 Mq (~ 100 Mq). Tables 5 and 6 report the frac¬ 
tion of BHs and massive stellar black holes (MSBHs, i.e. BHs 
with mass > 25 Mq, according to the definition by Mapelli 
et al. 2010) that form in our runs. 

The fraction of BHs in Table 5 is remarkably similar in 
all compared codes. Furthermore, this number is almost in¬ 
dependent of metallicity. On the other hand, the tested codes 
exhibit significant differences when the fraction of MSBHs 
is considered (Table 6). At metallicity Z = 2.0 x 10 2 , none 
of the compared codes form MSBHs, in agreement with the 
mass spectra we presented in Fig. 12. At lower metallic¬ 
ity, SEVN produces, on average, 5 — 6 times more MSBHs 
than SSE and StarlabMM. Therefore, the PARSEC stel¬ 
lar evolution prescriptions (combined with Fryer et al. 2012 
SN models) tend to form, approximately, the same number 
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Table 6. Fraction of MSBHs (that is, BHs with mass > 25 Mg), 
normalized to the total number of stars, obtained from SEVN, 
adopting different SN explosion models, and from StarlabMM 
and SSE. D: delayed model; R: rapid model; S: StarTrack pre¬ 
scriptions. 


z 


SEVN 


StarlabMM 

SSE 


D 

R . 

S 



2.0 x 10“ 4 

1.04 

1.00 

1.30 

0.20 

0.16 

2.0 x 10“ 3 

1.00 

0.96 

1.24 

0.18 

0 

2.0 x 10“ 2 

0 

0 

0 

0 

0 


The values are normalized to 10 3 . 


of BHs as the other codes, but many more MSBHs at low 
metallicity. 


4 COMPARISON WITH 

COMPACTNESS-BASED AND 
TWO-PARAMETER MODELS 

The Fryer et al. (2012) SN models we described in Sec¬ 
tions 3.3-3 .6 (as well as the other explosion prescriptions 
adopted in A-body simulations so far) are based on a 
single-parameter criterion that discriminates between SN 
explosion or failed SN. In this framework, stars explode if 
Mco < A/co,cut with A/co,cut = 11.0 M 0 for the rapid and 
delayed models and Mco,cut = 7.6 M 0 for the StarTrack 
model. Recent studies have shown that the link between 
physical properties of the progenitor star, SN properties 
and mass of the compact remnant is far from being trivial 
(Sukhbold & Woosley 2014; Ugliano et al. 2012; O’Connor 
& Ott 2011; Smartt 2015; Janka 2012; Ertl et al. 2015). 
In particular, it has been shown that the internal structure 
of the stars at core collapse varies non-monotonically with 
AFzams (or Mco) and this may lead to different compact 
remnants even if the progenitors had very similar A/zams- 
In this Section, we highlight the main differences between 
criteria based on Mco,cut and more sophisticated models, 
based on the structural properties of the star at the pre-SN 
stage. 

4.1 The compactness criterion 

O’Connor & Ott (2011) suggest that the value of the com¬ 
pactness £m evaluated just outside the iron core can discrim¬ 
inate between SNe and failed SNe. £m is the ratio between 
the innermost mass M of the star, in units of M 0 , and the 
radius R (M) containing M, in units of 1000 km, i.e. 

= M/M 0 

“ R(M )/1000 km' y ) 

Large values of £m favour failed SNe, while SNe occur for 
small values of £m- Generally, a fiducial value of M = 2.5 M 0 
is used to evaluate the compactness just outside the iron 
core. Even if the value of £ 2.5 is sensible to changes in mass 
loss prescriptions and stellar evolution parameters (such 
as mixing, reaction rates, opacity, metallicity), a threshold 
£ 2.5 ~ 0.2 seems to be a reasonable value to distinguish be¬ 
tween the occurrence of explosion and failed SNe (Horiuchi 
et al. 2011; Smartt 2015). Hereafter, we refer to the £ 2 . 5 - 
parameter model as £-modcl. 


Z = 2.0E-2, 2p-model 



Figure 23. Representation of the results we obtained with PAR¬ 
SEC progenitors in the two-parameter space introduced by Ertl 
et al. (2015), at Z = 0.02. Filled circles represent the formation 
of BHs via direct collapse, while open triangles refer to SNe. The 
dashed curve that divides SNe from failed SNe comes from cali¬ 
bration wl8.0 of Ertl et al. (2015). 


In order to use the £-model, SEVN needs further infor¬ 
mation (in addition to the standard input tables described 
in Appendix A), that is (i) the value of R(M) at the core 
collapse stage 10 , to evaluate £ 2.5 and to distinguish between 
SNe and failed SNe; (ii) the mass of the iron core Mp e , which 
is taken as the mass of the proto-compact object M pro to- 
Since PARSEC numerically integrates the stellar struc¬ 
ture up to the beginning of the CO burning phase only, we 
merged the PARSEC wind prescriptions with the MESA 
code (Paxton et al. 2011, 2013) and we used MESA to 
evolve the PARSEC models until the iron core infall phase. 
Our grid of MESA simulations goes from Mzams = 10 M 0 
up to A/zams = 30 M 0 with steps of 0.1 M 0 , and from 
A/zams = 30 M 0 to AFzams = 50 M 0 with steps of 0.3 M 0 . 

Fig. 21 shows how the compactness parameter £ 2.5 
changes as a function of Mzams- In this plot, we chose a 
critical compactness value of £2.5 = 0 .2 (Horiuchi et al. 2011, 
2014) to separate SNe from failed SNe. The relation between 
£ 2.5 and Mzams is quite complex. In particular, with the £- 
model, we distinguish at least three areas: 

(i) range Mzams £ [10 M 0 ; 18 M 0 ]: the majority of stars 
explode as SN and leave a NS with mass A/p e (excluding 
fallback material); 

(ii) range A/zams £ [18M 0 ;26M 0 ]: both SNe and failed 
SNe occur in this mass range; 

(iii) range A/zams > 26 M 0 : the majority of stars un¬ 
dergo direct collapse forming a BH with mass Mbh = Ma n . 

Even if the value of £ 2.5 depends on many stellar evolution 
parameters (including the adopted mass loss recipes), we 
find similar results to those obtained by other authors (e.g. 
Ugliano et al. 2012; Ertl et al. 2015). 


4.2 The two-parameter model 

A recent study by Ertl et al. (2015) introduces a two- 
parameter criterion. The two parameters are M4, which rep¬ 
resents the enclosed mass at a dimensionless entropy per nu¬ 
cleon s = 4, and fj. 4 , that is the mass gradient at the same 
location. Following the definition by Ertl et al. (2015), M 4 is 
normalized to M 0 and /n is normalized to 10 3 km/ M 0 . Ertl 
et al. (2015) show that a separation curve exists, that divides 


10 In our models, we identify the pre-SN stage when the collapse 
speed reaches ~ 10 s cm/s. 
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^ZAMS ( ) 


Figure 21. Value of compactness parameter at the innermost 2.5Mq, £ 2 . 5 ? as a function of Mzams? f° r PARSEC models evolved 
until the Fe-core infall using MESA. Black bars indicate non-exploding models (failed SNe) while white bars refer to exploding models 
(SNe). The dotted line £ 2.5 = 0-2 is the threshold we chose to distinguish between SNe and failed SNe according to Horiuchi et al. (2014). 
The simulation grid goes from Mzams — 10.0 Mq up to Mzams = 30.0 Mq with steps of 0.1 M©, and from Mzams = 30.0 Mq to 
Mzams = 50.0 Mq with steps of 0.3 Mq. Some models in the grid are not shown in the results because of numerical convergence issues. 



Figure 22. Same as Fig. 21 but here we show the parameter M 4 of the model by Ertl et al. (2015) as a function of Mzams- M 4 
represents the baryonic mass of the proto-compact object following the SN explosion event. In this case, to separate between SNe and 
failed SNe, we used the linear function y se p {%) = 0.283a: + 0.0430 which corresponds to the calibration wl8.0 of Ertl et al. (2015). 


exploding from non-exploding stars, in the plane x = M 4 /M, 
y = /X 4 . The threshold function is a straight line 

y se P (®) = fci x + fc 2 (8) 

where the coefficients k\ and fc 2 slightly depend on the differ¬ 
ent calibrations of the free parameters of Ertl et al. (2015) 
ID hydrodynamical simulations. Here, we use the calibra¬ 
tion curve for the model wl8.0 given by Ertl et al. (2015), 
for which fci = 0.283 and fc 2 = 0.043. Progenitors with 
J/pi-ogenitor > J/aep collapse directly into a BH, otherwise they 
explode as SN. Hereafter, we refer to this model as 2p-model. 
In order to apply this criterion to PARSEC progenitors, we 
extract the values of A /4 and fM from our grid of simulations 
run with MESA, coupled with the PARSEC wind models 
(see Sec. 4.1 for details). 

Fig. 22 shows the parameter A /4 (baryonic mass of the 
remnant if the compact object is a NS) as a function of 
A/zams ■ Black bars indicate direct collapse while white bars 
refer to SN explosion events. In this plot, we distinguish four 
different regions: 

(i) range A/zams £ [ 10 Mq; 18Mq]: the majority of stars 
explode as SNe and leave a NS with baryonic mass A /4 (ex¬ 
cluding fallback material); 

(ii) range A/zams € [18Mq; 24M 0 ]: both SNe and failed 
SNe occur; 

(iii) range A/zams £ [24 M 0 ; 28 M@]: the majority of stars 
undergo SN explosion; 

(iv) range A/zams > 28 M 0 : the majority of stars form 
BHs through direct collapse. 

For the calibration we assume, the main difference 


with the £-model (see Fig. 21) is in the range A/zams £ 
[24 M 0 ; 28 M 0 ], where the 2p-model produces a significantly 
higher number of NSs. This result confirms that BHs (NSs) 
can form even for A/zams < 25M 0 (A/zams > 25M 0 ). 

Fig. 23 shows the parameter y = im as a function of 
x = A/ 4 /M, for the 2p-model. Filled circles indicate BH for¬ 
mation (via direct collapse) while open triangles refer to the 
production of NSs (SNe explosion). Our PARSEC progeni¬ 
tors populate a narrow region in the x — y parameter space, 
whose range is similar to that shown in Ertl et al. (2015). 


4.3 Comparison with FFyer et al. (2012) models 

Fig. 24 shows the mass spectrum of compact remnants, ob¬ 
tained using the £-model (filled circles) and the 2p-model 
(open triangles), as a function of A/zams (black points), at 
Z = 0.02. In the same figure, we also represent the mass 
spectrum given by the delayed, rapid and StarTrack mod¬ 
els. Since both the £-model and the 2p-model do not provide 
prescriptions to evaluate the amount of mass that falls back 
onto the proto-compact object, all the models shown in Fig. 
24 do not include fallback. 

Overall, the mass spectrum of compact remnants re¬ 
sulting from either the £-model or the 2p-model is similar 
to the one derived from the StarTrack model. The main 
difference is that the £- and 2p-models predict a significant 
amount of BHs, due to failed SNe, for A/zams < 30 M 0 . 
Using the delayed and rapid models, direct collapse occurs 
for A/zams > 50 M 0 only. 

Finally, Fig. 24 shows that there is a significant mass 
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Figure 24. Mass of the compact remnant as a function of Mzams 
for the SN explosion recipes discussed in this paper, at Z = 0.02. 
In particular, filled circles: £—model; open triangles: 2p-model; 
solid line: delayed SN model; long-dashed line: rapid SN model; 
short-dashed line: Startrack recipes. In this plot, fallback is not 
included. We insert a y-axis break between 2 Mg and 11 Mg to 
better represent the mass spectrum of both BHs and NSs for the 
various models and the mass gap between the heaviest NS and 
the lightest BH. 

gap between the heaviest NS (~ 2 Mg) and the lightest BH 
(~ 12 Mg), quite larger than the observed one (Ugliano et al. 
2012; Farr et al. 2011; Ozel et al. 2010; Neustroev et al. 
2014; Kochanek 2014). Still, the area between 2 Mg and 12 
Mg may be populated by NSs that accrete mass through 
the fallback mechanism 11 , and/or by NSs that accrete mass 
from a companion, in a binary system. 


5 CONCLUSIONS 

The mass spectrum of BHs is still an open issue: only a few 
dynamical mass measurements of BHs are available (Ozel 
et al. 2010), while theoretical models are affected by the 
uncertainties on SN explosion and massive star evolution. 
In this paper, we derive the mass spectrum of compact 
remnants based on the new stellar evolution models im¬ 
plemented in PARSEC (Bressan et al. 2012; Chen et al. 
2014; Tang et al. 2014), combined with different recipes for 
SN explosion: the rapid and delayed SN models presented 
in Fryer et al. (2012), the SN model implemented in the 
StarTrack code (Belczynski et al. 2008), the SN recipes 
included in Starlab through the SeBa module (Portegies 
Zwart et al. 2001), and (for Z = 0.02) the £— and the 2p- 
models (O’Connor & Ott 2011; Ertl et al. 2015). 

These recipes for stellar evolution and SN explosion 
are implemented in our new public tool SEVN, which can 
be used both as a stand-alone population-synthesis code 

11 Ertl et al. (2015) show that a typical amount of mass that falls 
back onto the proto-compact object is ~ 0.05 Mg. Values larger 
than 1 Mgare rare (only 6 events over ~ 600 progenitor models). 


or as a module in several N-body codes (Starlab, HiG- 
PUs). SEVN is extremely versatile, because it calculates 
the mass, radius, luminosity, temperature and chemical evo¬ 
lution of a star based on stellar-evolution tables. We adopt 
stellar-evolution tables that have been generated with the 
PARSEC code, but these can be substituted with different 
stellar-evolution models in a fast and simple way. 

With respect to previous stellar-evolution codes, PAR¬ 
SEC predicts significantly larger values of Ma n and Mco 
at low metallicity (< 2 x 10~ 3 , Figures 4 and 5). We find 
differences up to ~ 80 % between the value of Ma n cal¬ 
culated by PARSEC and the fitting formulas implemented 
in SSE (Fig. 3). This implies that SEVN predicts substan¬ 
tially larger BH masses at low metallicity, since the mass of 
the compact remnants depends on Ma n and Mco in the SN 
models developed by Fryer et al. (2012). 

Moreover, for a metallicity Z = 0.02 and for Mzams 5) 
50Mg, we also present the mass spectrum of NSs and BHs 
given by the £-model (O’Connor & Ott 2011) and the 2p- 
model (Ertl et al. 2015). These models depend on stellar 
structural parameters evaluated at the time of iron core in¬ 
fall. Coupling these new prescriptions with the PARSEC 
stellar models, we find that the relation between progen¬ 
itor mass and remnant mass is quite complex, especially 
in the range Mzams £ [18 Mg; 30 Mg] (see Figures 21 and 
22). A detailed study that considers also Z ^ 0.02 and 
Mzams > 50 Mg is still in progress. 

Using the Fryer et al. (2012) models, we find that the 
maximum BH mass found with SEVN is ~ 25, 60 and 130 
Mg at Z = 2 x 10 2 , 2 x 10~ 3 and 2 x 10 —4 , respectively. 
Mass loss by stellar winds plays a major role in determining 
the mass of BHs for very massive stars (> 90 Mg), almost 
independently of the adopted SN recipe. In contrast, the 
adopted SN model is very important for lower BH masses, 
and for the transition between NSs and BHs (Figures 8, 9, 
10 and 11): according to the delayed SN model, stars with 
Mzams > 19 Mg end their life as BHs, while this limit is 
Mzams > 24 — 25 Mg if the rapid SN mechanism or the 
SeBa recipes are assumed. 

As a consequence, the rapid SN mechanism and the 
recipes implemented in SeBa predict a gap between the 
maximum mass of NSs and the minimum mass of BHs, while 
the delayed SN model (and the StarTrack recipes) suggest 
a smooth transition between NSs and BHs (Figs. 15, 16 and 
17). The distribution of dynamically measured BH and NS 
masses in the local Universe suggests the existence of a gap 
between NS and BH masses (Ozel et al. 2010), even if the 
statistical significance of this result is still debated (Farr 
et al. 2011; Ugliano et al. 2012; Kochanek 2014; Neustroev 
et al. 2014). 

According to SEVN (with either the delayed or the 
rapid SN model), at Z = 2 x 10 -2 most BHs have mass 
8 — 12 Mg, while at 2 x 10~ 3 ^ 2 ^ 2 x 10 -4 most BHs 
have mass 20 — 60 Mg (Figs. 15, 16 and 17). 

For a stellar population following the Kroupa IMF, the 
total number of BHs predicted by SEVN in its various SN 
flavours is remarkably similar to other codes, such as Star- 
labMM (Mapelli et al. 2013) and SSE (Hurley et al. 2000). 
Furthermore, the fraction of BHs is almost independent of 
metallicity. On the other hand, the fraction of MSBHs (i.e. 
BHs with mass > 25 Mg) strongly depends on the metallic¬ 
ity and on the assumed stellar evolution recipes. At metallic- 
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ity Z = 2.0 x 10 2 , no MSBHs form from single-star evolu¬ 
tion, using either SEVN or StarlabMM or SSE. At lower 
metallicity, SEVN produces, on average, 5 — 6 times more 
MSBHs than SSE and StarlabMM. 

This might have dramatic consequences for both the 
number of X-ray binaries powered by MSBHs and the de¬ 
tection of gravitational waves by BH-BH binary mergers. 
As to X-ray binaries, models by Mapelli & Zampieri (2014), 
based on StarlabMM, indicate that MSBHs are expected 
to power ~ 20 % of the Roche-lobe overflow BH binaries in 
a young star cluster with Z < 2 x 10 -3 . With the recipes 
implemented in SEVN, the fraction of X-ray binaries pow¬ 
ered by MSBHs might be substantially higher. On the other 
hand, quantifying the difference with previous studies is non 
trivial, because the evolution of binary systems and dynam¬ 
ical encounters in star clusters can significantly affect the 
demographics of BH binaries. In a forthcoming study, we 
will use SEVN to investigate the demographics of X-ray bi¬ 
naries and BH-BH binaries in star clusters. 
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APPENDIX A: NOTES ON THE 
IMPLEMENTATION OF SEVN 

A1 General scheme 

As we discussed in Sec. 2.2, SEVN can be used as a stand¬ 
alone population synthesis code and/or can be easily linked 
to several TV-Body codes. It is extremely versatile because 
it relies upon a set of isochrones as input files; this means 
that we can easily change stellar evolution recipes by simply 
substituting the input tables. 

SEVN reads a single file that contains a set of 
isochrones, which are provided by a stellar evolution code 
(by PARSEC, in the current version of SEVN). By default, 
the name of this file must have the form tableZn.dat, where 
n indicates the nretallicity Z. Inside this file, the isochrones 
are separated by a line reporting the age (in Gyr), and the 
number of points of the following isochrone. Each isochrone 
is composed of 9 columns that indicate (1) the initial mass 
of the star, (2) its present mass, (3) the logarithm of lumi¬ 
nosity, (4) the effective temperature, (5) the logarithm of 
radius and (6) the logarithm of surface gravity, (7) the He¬ 
lium and (8) Carbon-Oxygen core mass and (9) the stellar 
type at the current age. All given values are in solar units, 
except for the effective temperature, which is absolute and 
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expressed in Kelvin. We stress that the isochrones do not 
need to be equally spaced in mass or other quantities. 

In order to speed up the calculations, SEVN reads the 
isochrone file and rearranges it in a more convenient way. 
First of all, an equally spaced grid of masses is chosen 12 . 
For each star in the grid, we construct the time evolution 
of its physical parameters, recording information whenever 
the value of a generic stellar parameter is varied by more 
then 5%. The result is stored in 7 different files contain¬ 
ing the time evolution of masses, radii, luminosities, stellar 
phases, Carbon-Oxygen core mass, Helium core mass and 
the corresponding ages when the stellar parameters need to 
be updated. 

These 7 files are then loaded in a 3-dimensional struc¬ 
ture where the first index (line number, L) identifies the 
initial mass of the star. The second index (column number, 
C) gives information about the current stellar age and the 
third index, P, refers to the specific stellar parameter we 
need to read or write. Thus, L ranges between 1 and the 
number of points of the grid of masses, 1 ^ P ^ 7, and C 
varies from 1 to the number of update points needed for a 
generic star. 

At the beginning of the integration, it is possible to 
associate two different mass indexes, L\ and L 2 , to each 
star in order to uniquely identify its position in the grid. 
For example, let us consider a grid of masses that goes from 
0.1 Mg to 150 Mg with steps of 0.5 Mg. The evolution of a 
star S of mass M s = 50.3MQwill be derived interpolating 
the evolutionary tracks of the nearest neighbour stars, that 
is Mi = 5OM0and M 2 = 50.5 Mg, and we can compute the 
stellar parameters of the star S using the weights 


M 2 - M a 

(Ms - Mi) + (M 2 - Ms) 
Ms - Mi 

(Ms - Mi) + (Ma -M s )' 


(Al) 


To evolve the parameters of a generic star, we use linear 
interpolations. Let us consider again a test star S with initial 
mass M s (t = 0) = 50.3 Mg. At time t = t\ , this star will 
have a mass M s (ti). In order to evolve the star at time 
t 2 = ti + Al, we need to use the information of its neighbour 
grid stars of mass Mi (0) = 50Mgand M 2 (0) = 50.5 Mg. 
First of all, the code must compute the quantities Mi (t 2 ) 
and M 2 (f 2 ). In general, a generic time f 2 will not be included 
in the tables. Thus, SEVN reads the tables and searches 
the values Mi ( 13 ), Mi ( 14 ), M2 (Is) and M2 (te) such that 
£3 ;$ t2 < ti and ts < t 2 < te- The code then calculates 
Mi (t 2 ) and M 2 (t 2 ) with a linear interpolation: 


Mi (t 2 ) = mif 2 + qi 
M 2 (1 2 ) = m 2 f 2 + q 2 


(A2) 


12 By default, the grid goes from 0.1 Mg to 150 Mg with steps 
of 0.5 Mg. 


where 


Mi (14) — Mi (tz) 

ti — I3 

M2 (te) — M2 (is) 
te — te 
qi = Mi (U) - miti 
q 2 = M 2 (te) - m,2te- 


m 1 = 


m 2 = 


(A3) 


Finally, the value M s (t 2 ) is derived with a further linear 
interpolation, that is 


Ms (ti) = a\M\ (t 2 ) + a 2 AL 2 (t 2 ) (A4) 

with weights on and a 2 given in equation Al. The same 
procedure is adopted to obtain the other stellar parameters 
needed at a given age. 


A2 Integration of SEVN in Starlab 

As discussed in Section 2.2, we have merged SEVN with 
the Starlab software environment (Portegies Zwart et al. 
2001), and with the direct V-Body code HlGPUs (Capuzzo- 
Dolcetta et al. (2013); Spera, in preparation). In order to 
combine SEVN with Starlab, we modified the SeBa stel¬ 
lar evolution module. In particular, SeBa is a C++ module 
based on a structure of classes, in which each class approx¬ 
imately corresponds to a stellar evolution phase. In order 
to easily match the SeBa internal organization and the im¬ 
plemented transitions between stellar evolution phases, we 
identify the main stellar evolution phases using integer in¬ 
dexes. Namely, 

• 0 identifies pre-MS and MS stars that are mapped to 
the SeBa class main_sequence; 

• 1 indicates stars in the sub-giant phase; in this case, we 
have a one to one correspondence with the sub_giant class; 

• 2 groups several categories of stars (among which, red 
giants, blue and red super giants, LBVs and WRs) in the 
hyper_giant class; 

• 3 refers to core helium burning stars, collected in the 
horizontal_branch class; 

• 4 corresponds to stars in the early asymptotic giant 
branch (E-AGB) phase, mapped to the hyper_giant class. 

Using our simplified scheme, we lose information about 
some specific characteristics of the stars during the numer¬ 
ical integration; for instance, we do not know if a star is 
a WR, a LBV or a blue super giant or a red super giant. 
Anyway, all these features can be recovered a posteriori by 
the ages, radii, luminosities and temperatures printed in the 
output files. 

During the pre-MS and MS phases we evolve mass, lu¬ 
minosity and radius of the stars following our input tables 
by means of linear interpolations in time and mass. Stel¬ 
lar evolution continues until the function create_remnant 0 
is called. This routine contains our updated recipes for SN 
explosion, and converts the star into a compact remnant, 
which can be either a WD, a NS or a BH, depending on the 
final state of the star. 

At present, since PARSEC does not include evolu¬ 
tionary prescriptions for stars that undergo the thermally- 
pulsing AGB phase (TP-AGB), we use the SeBa built-in 
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class super_giant to follow their evolution through this 
stage (see Portegies Zwart et al. (2001) for the details). In 
particular, we assume that the stars that undergo the TP- 
AGB phase are those with Mzams < M up = 7Mq. 


APPENDIX B: SN EXPLOSION MECHANISMS 
IN SEVN 

Here we summarize the main features of the Fryer et al. 
(2012) recipes. 


BO. 1 StarTrack model 

In the case of StarTrack recipes, stars form a proto¬ 
compact object of mass M proto given by 


0.2 

-Mfi n Al proto 

0.286Mco - 0.514 

A/n n A/j lro to 

1.0 

orMco + /3r 

,1.0 


where 


Mco < 2.5 Mg 


2.5 Mg ^ Mco < 6.0 Mg 


6.0 Mg ^ Mco < 7.0 Mg 
7.0 Mg ^ Mco < 11.0 M© 
Mco ^ 11.0 M© 

(B5) 


= 0.25 — 


1.275 

ALfl n Mproto 


Bit = 1 - lla R . 


(B6) 


This means that the direct collapse of a star into a BH 
occurs if 6.0 M© ^ Mco ^ 7.0 M© and if Mco 11 M© 
(equation B5), according to the rapid SN model. 


Mr 


proto : 


' 1.50 M© 

2.11 M© 

| 0.69AL C o -2.26M© 
,0.37AL C o -0.07M© 


Mco < 4.82 Mq 
4.82 Mq ^ Mco < 6.31 M© 
6.31 M© < Mco < 6.75 M 0 
Mco ^ 6.75 Mq. 


(Bl) 

/fb is the fractional fallback parameter , and is such that 
Mr, = /fb (Mfl n — Mp roto ), where Affl n is the final mass of 
the star. According to StarTrack prescriptions, the values 
of /fb are the following: 


[0 Mco < 5.0 M© 

/fb = < 0.378Af C o - 1-889 5.0 M© sC M C o < 7.6 M© 
[l.O Mco ^ 7.6 M©. 

(B2) 

From the baryonic mass of the remnant M rern ^ ar = M pr oto + 
ALfb, we can obtain its gravitational mass Mrem,gray taking 
into account neutrino losses. When Mco ^ 7.6 M©, the 
StarTrack recipes request /n, = 1 (in eq. B2), i.e. the 
entire final mass of the star goes into the remnant mass. This 
means that the direct collapse of a star into a BH occurs if 
Mco ^ 7.6 M©, according to StarTrack prescriptions. 

For NSs we use the expression given by Timmes et al. 
(1996), for which 


BO. 3 Delayed SN model 

For the delayed SN mechanism, the prescriptions for the 
mass of the proto-compact object are 


' 1.2 M© Mco < 3.5 M© 

1.3 M© 3.5 M© ^ Mco < 6.0 M© 

1.4 M© 6.0 M© ^ Mco < 11.0 M© 
, 1.6 M© Mco ^ 11.0 M©. 


(B7) 


The amount of fallback is determined using the following 
relations 


fib 


0.2 

M fl,, Mp r0 t 0 

0.5M CO - 1.05 M© 
' M fln - Mp ro to 
CUD Mco + 0D 
1.0 


where 


Mco < 2.5 M© 

2.5 M© ^ Mco < 3.5 M© 

3.5 M© ^ Mco < 11.0 M© 
Mco ^ 11.0 M© 

(B8) 


CtD 


0.133 


0.093 

A f fi n Alproto 


0D = 1 — 11 OLD- 


(B9) 


Mrem.grav 


\/1 + 0.3A/,.,, rn ,bar 1 

0T5 


For BHs we use the formula 


(B.3) 


Thus, the direct collapse of a star into a BH occurs if 
Mco M© (equation B8), according to the delayed SN 
model (i.e. the same as the rapid SN model, but significantly 
larger than in the StarTrack recipes). 


Ml re m,grav — 0.9A4 rem ,6ar, (B4) 

following the approach described in Fryer et al. (2012). 

BO. 2 Rapid SN model 

For the rapid SN mechanism, a fixed mass of the proto¬ 
compact object, Mp rci to = 1.0 M©, is assumed. In this case, 
the coefficient /fb is given by 


APPENDIX C: GENERAL FITTING FORMULA 

FOR M rem 

We report a fitting formula that express the compact rem¬ 
nant mass Mrem as a function of Mzams and Z. The fol¬ 
lowing formula has been obtained by fitting the outputs of 
SEVN with the delayed SN model and the PARSEC stel¬ 
lar evolution isochrones. The value of M rem obtained with 
the fitting formula deviate from the outputs of SEVN by 
< 10 %. 
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First, we express A/ rem as a function of Mco and Z 
(from Fig. 7). For Z ^ 5.0 x 10~ 4 , the best fitting curve for 
Mrem is given by: 


('max (p (Mco), 1.27M 0 ) 
if Mco ^ 5 M 0 


M rem = < 


P (Mco) 

if 5 Mg < Mco < 10 Mg 


(Cl) 


min (p ( Mco ), / (M C o, Z)) 
k if Mco ^ 10Mg, 

where 

P (M CO ) = -2.333 + 0.1559M CO + 0.2700A/co 
f (Mco, Z) = m (Z) Mco + q (Z) 
with coefficients 

m (Z) = -6.476 x 10 2 Z + 1.911 
q(Z) = 2.300 x 10 3 Z+ 11.67. 

For Z > 5.0 x 10~ 4 we have 

max (h (Mco, Z), 1.27Mg) 
if Mco ^ 5 M 0 

h(M CO ,Z) 

if 5 Mg < Mco ^ 10 M 0 

max (h (Mco, Z), f (M C o,Z)) 
if Mco ^ 10 Mg, 


(C2) 


(C3) 


M rem = < 


(04) 


where 


M 2 (Z) - Ai (Z) 


h (Mco, Z) Mi (Z) + 1 + 10 (Ii(Z )_m C o)„(z) 
f (Mco, Z) = m (Z) Mco + 9 (Z). 


(05) 


For Z ^ l.Ox 10 3 , the coefficients of the function 
h (Mco, Z) are 

Mi (Z) = 1.340 - 29 ' 46 


1 + 

M 2 (Z) = 80.22 - 74.73 
L (Z) = 5.683 + — 


1.110 x lO" 3 

^0.965 


2.720 x 10-3 + Z 0 - 965 
3.533 


(C6) 


1 + 


p(Z) = 1.066 - 


7.430 x lO" 3 
1.121 


1 + 


2.558 x 10- 


while for Z < 1.0 x 10 , the coefficients of the function 

h (Mco, Z) are 


Mi (Z) = 1.105 x 10 5 Z - 1.258 x 10 2 
M 2 (Z) = 91.56 - 1.957 x 10 4 Z - 1.558 x 10 7 Z 2 
L (Z) = 1.134 x 10 4 Z - 2.143 
p (Z) = 3.090 x 10“ 2 - 22.30Z + 7.363 x 10 4 Z 2 . 


( 07 ) 


For Z Js 2.0 x 10 3 , the coefficients of the function 
/ (Mco, Z) are independent of Z, 


m = 1.217 
q = 1.061, 

while, for 1.0 x 10~ 3 ^ Z < 2.0 x 10 -3 we have 


(C8) 


m= -43.82Z+ 1.304 
q = -1.296 x 10 4 Z + 26.98, 
and for Z < 1.0 x 10 -3 


(C9) 


m = -6.476 x 10 2 Z+ 1.911 
q = 2.300 x 10 3 Z+ 11.67. 


(CIO) 


Furthermore, A/co can be expressed as a function of 
A/zams and Z, by fitting the curves of Fig. 5. The functional 
form of the fit is 


Mco = —2.0 + [Si ( Z ) + 2.0] [g (Z, A/zams; K\, 8 i) + 


+ g (Z, A/zams; K2 , <5 2 )], 


(Cll) 


where 


g (Z, Mzams; x, y) - 1 + W (x(z)-M ZAMS ) y (z) ■ ( C12 ) 

For Z > 4.0 x 10 -3 the coefficients are 

Si (Z) = 59.63 - 2.969 x 10 3 Z + 4.988 x 10 4 Z 2 
AT (Z) = 45.04 - 2.176 x 10 3 Z + 3.806 x 10 4 Z 2 
A 2 (Z) = 1.389 x 10 2 - 4.664 x 10 3 Z + 5.106 x 10 4 Z 2 
<5i (Z) = 2.790 x 10" 2 - 1.780 x 10 _2 Z + 77.05Z 2 
82 (Z) = 6.730 x 10" 3 + 2.690Z - 52.39Z 2 . 

(C13) 

For 1.0 x 10 -3 ^ Z ^ 4.0 x 10 -3 , we have 

Si (Z) = 40.98 + 3.415 x 10 4 Z - 8.064 x 10 6 Z 2 
AT (Z) = 35.17 + 1.548 x 10 4 Z - 3.759 x 10 6 Z 2 
AT (Z) = 20.36 + 1.162 x 10 5 Z - 2.276 x 10 7 Z 2 (C14) 

Ai (Z) = 2.500 x 10 -2 - 4.346Z + 1.340 x 10 3 Z 2 
82 (Z) = 1.750 x 10 -2 + 11.39Z - 2.902 x 10 3 Z 2 . 

Finally, for Z < 1.0 x 10~ 3 , the coefficients do not de¬ 
pend on Z, 

Si = 67.07 
AT =46.89 

AT = 1.138 x 10 2 (C15) 

Ai = 2.199 x 10~ 2 
82 = 2.602 x 10~ 2 . 



